This notebook draws from previous notebooks on relative cumulative frequency analysis and PCA analysis of metabolic cage data.
These datasets are from RMR data collected for mice after introduction to the metabolic cages (days 0-6, see image).
A total of 24 mice were run over the length of 2 experiments.
Experiment 1: The first set of mice, run in February
| Start of Met Cage run |
Metabolic Cage # |
Ear Tag | DOB | Gender | Genotype | Body Mass on BMR day (g) |
|---|---|---|---|---|---|---|
| 2/4/2020 | 1 | 1049 | 10/08/19 | M | WT | 27.4 |
| 2/4/2020 | 2 | 1050 | 10/08/19 | M | L-Bmal1-KO | 26.6 |
| 2/4/2020 | 3 | 1053 | 10/07/19 | M | WT | 31.5 |
| 2/4/2020 | 4 | 1054 | 10/07/19 | M | WT | 30.5 |
| 2/4/2020 | 5 | 1055 | 10/07/19 | M | L-Bmal1-KO | 27.5 |
| 2/4/2020 | 6 | 1109 | 10/11/19 | M | WT | 28.3 |
| 2/4/2020 | 9 | 1110 | 10/11/19 | M | WT | 27.9 |
| 2/4/2020 | 10 | 1111 | 10/11/19 | M | L-Bmal1-KO | 26.7 |
| 2/4/2020 | 11 | 1112 | 10/11/19 | M | WT | 28.5 |
| 2/4/2020 | 12 | 1214 | 10/21/19 | M | L-Bmal1-KO | 25.3 |
| 2/4/2020 | 13 | 1215 | 10/21/19 | M | WT | 27.9 |
| 2/4/2020 | 14 | 1216 | 10/21/19 | M | L-Bmal1-KO | 27.6 |
Experiment 2: The second set of mice, run in March
| Start of Met Cage run |
Metabolic Cage # |
Ear Tag | DOB | Gender | Genotype | Body Mass on BMR day (g) |
|---|---|---|---|---|---|---|
| 3/10/2020 | 1 | 1555 | 11/18/19 | M | WT | 26.4 |
| 3/10/2020 | 2 | 1692 | 12/06/19 | M | WT | 25.6 |
| 3/10/2020 | 3 | 1693 | 12/06/19 | M | WT | 24.7 |
| 3/10/2020 | 4 | 1695 | 12/06/19 | M | WT | 25.6 |
| 3/10/2020 | 5 | 1699 | 12/06/19 | M | WT | 28.7 |
| 3/10/2020 | 6 | 1556 | 11/18/19 | M | L-Bmal1-KO | 26.2 |
| 3/10/2020 | 9 | 1557 | 11/18/19 | M | L-Bmal1-KO | 27.1 |
| 3/10/2020 | 10 | 1558 | 11/18/19 | M | L-Bmal1-KO | 26.9 |
| 3/10/2020 | 11 | 1559 | 11/18/19 | M | L-Bmal1-KO | 24.7 |
| 3/10/2020 | 12 | 1691 | 12/06/19 | M | L-Bmal1-KO | 28 |
| 3/10/2020 | 13 | 1694 | 12/06/19 | M | L-Bmal1-KO | 26.4 |
| 3/10/2020 | 14 | 1700 | 12/06/19 | M | L-Bmal1-KO | 27.7 |
Data generated using a python script.
SPF_vs_GF.py
#!/usr/bin/env python
"""
A CLI python script that uses the mousetrapper library to prepare WT vs KO mice data for import to R.
Must be run from same directory.
"""
import argparse
import numpy as np
import pandas as pd
import lib.mousetrapper as mt
__author__ = "Sumeed Yoyo Manzoor"
__email__ = "smanzoor@uchicago.edu"
__version__ = "0.0.1"
__experiment__ = "12 WT vs 12 L-Bmal1-KO mice"
datasets = "../datasets/"
data = "../data/"
outfolder = data + "R/"
experiment1_RMR = mt.Experiment.experiment(
info="""First run with revised protocol, WT vs L-Bmal1-KO""",
RT_data=datasets + "Bmal_WT_vs_KO_SPF_RMR_2020-02-04_rt.csv",
macro13_data=datasets + "bmal_wt_vs_ko_spf_rmr_2020-02-04_macro_One-Click Macro V2.32.2-slice3-VOC.mac_1.csv",
t_ctrl={"SPF WT":[1,3,4,6,9,11,13]},
t_exp={"SPF L-Bmal1-KO":[2,5,10,12,14]},
groups={
"SPF WT":[1,3,4,6,9,11,13],
"SPF L-Bmal1-KO":[2,5,10,12,14]
},
meta=["WT","L-Bmal1-KO","WT","WT","L-Bmal1-KO","WT","WT","L-Bmal1-KO","WT","L-Bmal1-KO","WT","L-Bmal1-KO"],
covar=[27.4,26.6,31.5,30.5,27.5,28.3,27.9,26.7,28.5,25.3,27.9,27.6]
)
experiment2_RMR = mt.Experiment.experiment(
info="""Second run with revised protocol, WT vs L-Bmal1-KO""",
RT_data=datasets + "Bmal_SPF_KF_RMR_2020-03-10_14_23_rt.csv",
macro13_data=datasets + "bmal_spf_kf_rmr_2020-03-10_14_23_macro_One-Click Macro V2.32.2-slice3-VOC.mac_1.csv",
t_ctrl={"SPF WT":np.arange(1,6).tolist()},
t_exp={"SPF L-Bmal1-KO":[6] + np.arange(9,15).tolist()},
groups={
"SPF WT":np.arange(1,6).tolist(),
"SPF L-Bmal1-KO":[6] + np.arange(9,15).tolist()
},
meta=np.array(["WT","L-Bmal1-KO"]).repeat([5,7]).tolist(),
covar=[26.4,25.6,24.7,25.6,28.7,26.2,27.1,26.9,24.7,28,26.4,27.7]
)
experiment1_RMR.macro13_lightcycles(days=2)
experiment2_RMR.macro13_lightcycles(days=2)
experiment1_RMR.macro13_dark.to_csv(outfolder + "1_dark.csv", index=False)
experiment1_RMR.macro13_light.to_csv(outfolder + "1_light.csv", index=False)
pd.DataFrame({"covar":experiment1_RMR.covar}).to_csv(outfolder + "1_covar.csv", index=False)
experiment2_RMR.macro13_dark.to_csv(outfolder + "2_dark.csv", index=False)
experiment2_RMR.macro13_light.to_csv(outfolder + "2_light.csv", index=False)
pd.DataFrame({"covar":experiment2_RMR.covar}).to_csv(outfolder + "2_covar.csv", index=False)
shell("conda activate met_cages_development && cd py && python Bmal.py")
reqpkg <- c("magrittr","dplyr","ggplot2","ggpubr","DT","reshape2")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "magrittr"
## [1] '1.5'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
## [1] "DT"
## [1] '0.13'
## [1] "reshape2"
## [1] '1.4.3'
theme_set(theme_pubr())
select <- dplyr::select
path <- "data/R/"
exp1 <- read.csv2(paste0(path, "1_total.csv"), sep = ",") %>%
set_colnames(paste0("exp1_",colnames(.)))
exp1_dark <- read.csv2(paste0(path, "1_dark.csv"), sep = ",") %>%
set_colnames(paste0("exp1_",colnames(.))) %>% slice(-c(241, 242, 483, 484))
exp1_light <- read.csv2(paste0(path, "1_light.csv"), sep = ",") %>%
set_colnames(paste0("exp1_",colnames(.)))
exp1_covar <- read.csv2(paste0(path, "1_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
exp2 <- read.csv2(paste0(path, "2_total.csv"), sep = ",") %>%
set_colnames(paste0("exp2_",colnames(.)))
exp2_dark <- read.csv2(paste0(path, "2_dark.csv"), sep = ",") %>%
set_colnames(paste0("exp2_",colnames(.))) %>% slice(-c(241, 482))
exp2_light <- read.csv2(paste0(path, "2_light.csv"), sep = ",") %>%
set_colnames(paste0("exp2_",colnames(.))) %>% slice(-c(241, 482))
exp2_covar <- read.csv2(paste0(path, "2_covar.csv"), sep = ",") %>%
mutate_all(function(x) as.numeric(as.character(x)))
covar <- c(exp1_covar[c(1,3,4,6,7,9,11),], exp2_covar[1:5,], exp1_covar[c(2,5,8,10,12),], exp2_covar[6:12,])
total2 <- exp1 %>%
bind_cols(exp2) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
light <- exp1_light %>%
bind_cols(exp2_light) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
dark <- exp1_dark %>%
bind_cols(exp2_dark) %>%
select(-contains(c("Enviro", "Time"))) %>%
mutate_all(function(x) as.numeric(as.character(x)))
total <- bind_rows(light, dark)
Select a tab above to view
light %>% select(contains("VO2")) %>% datatable(filter="bottom")
light %>% select(contains("RQ")) %>% datatable(filter="bottom")
light %>% select(contains("BodyMass")) %>% datatable(filter="bottom")
light %>% select(contains("Food")) %>% datatable(filter="bottom")
light %>% select(contains("Water")) %>% datatable(filter="bottom")
light %>% select(contains("PedMeters")) %>% datatable(filter="bottom")
Select a tab above to view
dark %>% select(contains("VO2")) %>% datatable(filter="bottom")
dark %>% select(contains("RQ")) %>% datatable(filter="bottom")
dark %>% select(contains("BodyMass")) %>% datatable(filter="bottom")
dark %>% select(contains("Food")) %>% datatable(filter="bottom")
dark %>% select(contains("Water")) %>% datatable(filter="bottom")
dark %>% select(contains("PedMeters")) %>% datatable(filter="bottom")
select_WT <- function(df) {
d1 <- df %>% select(contains("exp1"))
for (name in colnames(d1)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(2,5,10,12,14)) {
d1 %<>% select(-all_of(name))
}
}
d2 <- df %>% select(contains("exp2"))
for (name in colnames(d2)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) > 5) {
d2 %<>% select(-all_of(name))
}
}
d <- bind_cols(d1, d2)
return(d)
}
select_KO <- function(df) {
d1 <- df %>% select(contains("exp1"))
for (name in colnames(d1)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,3,4,6,9,11,13)) {
d1 %<>% select(-all_of(name))
}
}
d2 <- df %>% select(contains("exp2"))
for (name in colnames(d2)) {
if (as.numeric(strsplit(name, "_")[[1]][4]) <= 5) {
d2 %<>% select(-all_of(name))
}
}
d <- bind_cols(d1, d2)
return(d)
}
graph_groups <- function(data, channel, title="", type="m") {
df <- data %>% select(contains(channel))
df.observation <- 1:nrow(df)
c <- ncol(df)
WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt() %>% select(-variable)
KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt() %>% select(-variable)
df <- bind_cols(WT, KO) %>% set_colnames(c("WT","KO")) %>% melt()
df$observation <- rep(df.observation, c)
if (type=="m") {
df %>% group_by(variable, observation) %>%
summarise(value=mean(value)) %>%
ggplot(aes(x=observation, y=value, color=variable)) +
geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_line(size=1) +
ggtitle(title) +
ylab(channel)
} else if (type=="s") {
ggplot(df, aes(x=observation, y=value, color=variable)) +
ggtitle(title) +
geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", alpha = 0.7) +
ylab(channel)
}
}
graph_mice <- function(data, channel, title="") {
df <- data %>% select(contains(channel))
df.observation <- 1:nrow(df)
c <- ncol(df)
WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt()
KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt()
df <- bind_rows(WT, KO) %>% set_colnames(c("mouse", "channel"))
df$observation <- rep(df.observation, c)
ggplot(df, aes(x=observation, y=channel, color=mouse)) +
geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
geom_line() +
ylab(channel)
}
Select a tab above to view
graph_groups(total2, "VO2", title="VO2 by group", type = "m")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_groups(total2, "VO2", title="VO2 by group", type = "s")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_mice(total2, "VO2", title="VO2 by mouse")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
graph_groups(total2, "RQ", title="RQ by group", type = "m")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_groups(total2, "RQ", title="RQ by group", type = "s")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_mice(total2, "RQ", title="RQ by mouse")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
graph_groups(total2, "BodyMass", title="BodyMass by group", type = "m")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_groups(total2, "BodyMass", title="BodyMass by group", type = "s")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_mice(total2, "BodyMass", title="BodyMass by mouse")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
graph_groups(total2, "Food", title="Food by group", type = "m")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_groups(total2, "Food", title="Food by group", type = "s")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_mice(total2, "Food", title="Food by mouse")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
graph_groups(total2, "Water", title="Water by group", type = "m")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_groups(total2, "Water", title="Water by group", type = "s")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_mice(total2, "Water", title="Water by mouse")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
graph_groups(total2, "PedMeters", title="PedMeters by group", type = "m")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_groups(total2, "PedMeters", title="PedMeters by group", type = "s")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables
graph_mice(total2, "PedMeters", title="PedMeters by mouse")
## No id variables; using all as measure variables
## No id variables; using all as measure variables
## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector
reqpkg <- c("docstring","tibble","scales","purrr","nplr", "car","compute.es","effects","multcomp","pastecs","WRS2")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "docstring"
## [1] '1.0.0'
## [1] "tibble"
## [1] '2.1.3'
## [1] "scales"
## [1] '1.1.0'
## [1] "purrr"
## [1] '0.3.3'
## [1] "nplr"
## [1] '0.1.7'
## [1] "car"
## [1] '3.0.7'
## [1] "compute.es"
## [1] '0.2.4'
## [1] "effects"
## [1] '4.1.4'
## [1] "multcomp"
## [1] '1.4.12'
## [1] "pastecs"
## [1] '1.3.21'
## [1] "WRS2"
## [1] '1.0.0'
signif.floor <- function(x, n){
pow <- floor( log10( abs(x) ) ) + 1 - n
y <- floor(x / 10 ^ pow) * 10^pow
# handle the x = 0 case
y[x==0] <- 0
y
}
signif.ceiling <- function(x, n){
pow <- floor( log10( abs(x) ) ) + 1 - n
y <- ceiling(x / 10 ^ pow) * 10^pow
# handle the x = 0 case
y[x==0] <- 0
y
}
find_range <- function(data, channel, segments=10) {
data <- data %>% dplyr::select(contains(channel)) %>% melt()
r <- range(data$value)
cut <- max(data$value)-min(data$value)
out <- list("range"=r, "cut"=cut/segments)
return(out)
}
rcf <- function(data, channel, id, min, max, b=0.1) {
#' Relative cumulative frequency function
#'
#' Takes a vector and returns a dataframe with relative cumulative frequency per break unit
#'
#' @param data numeric vector
#' @param channel character. Channel to label the rcf output.
#' @param b double. Number to cut vector.
#' Before running, check for valid break values to cut the data. One way is to divide range of data by 10 `(max(data)-min(data))/10`.
breaks <- seq(min, max, by=b)
# breaks <- seq(signif.floor(min(data), 2), signif.ceiling(max(data), 2), by=b)
duration.cut <- cut(data, breaks, right = F)
duration.freq <- table(duration.cut)
duration.cumfreq = cumsum(duration.freq) %>% prepend(0)
duration.cumrelfreq = duration.cumfreq/length(data)
out <- duration.cumrelfreq %>% as.data.frame()
out <- cbind(breaks, out) %>% set_colnames(c("breaks", channel)) %>% set_rownames(NULL)
out$id <- rep(id, nrow(out))
return(out)
}
EC50 <- function(x) {
out <- lapply(x, function(y) {
nplr(y$breaks, y[, 2], silent = T, useLog = FALSE)
}) %>%
sapply(function(z) {
getEstimates(z, 0.5)$x
})
out %<>% set_rownames(NULL)
return(out)
}
lm_eqn <- function(df, y, x){
m <- lm(y ~ x, df);
eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,
list(a = format(unname(coef(m)[1]), digits = 2),
b = format(unname(coef(m)[2]), digits = 2),
r2 = format(summary(m)$r.squared, digits = 3)))
as.character(as.expression(eq));
}
lm_eqn_str <- function(df, y, x){
m <- lm(y ~ x, df);
eq <- sprintf("y = %s + %sx, r^2 = %s", format(unname(coef(m)[1]), digits = 2),format(unname(coef(m)[2]), digits = 2),format(summary(m)$r.squared, digits = 3))
return(eq)
}
remove_outliers <- function(x, na.rm = TRUE, ...) {
qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
H <- 1.5 * IQR(x, na.rm = na.rm)
y <- x
y[x < (qnt[1] - H)] <- NA
y[x > (qnt[2] + H)] <- NA
y
}
ANCOVA <- function(df) {
# Takes dataframe with column 1 as independent variable (e.g., WT), column 2 as dependent variable, and column 3 as covariate
cat("linear fit of dependent variable with covariates\n")
p1 <- ggplot(df, aes_string(x=names(df)[2], y=names(df)[3])) +
geom_smooth(method = "lm", se = FALSE) +
geom_point() +
# geom_text(x = mean(na.omit(df[[2]])), y = mean(na.omit(df[[3]])), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
geom_text(x = mean(na.omit(df[[2]])), y = (mean(na.omit(df[[3]])) + IQR(df[[3]], na.rm = T)), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
ggtitle("Linear model - DV and Covar") +
xlab("DV") + ylab("Covar")
print(p1)
m <- lm(df[[3]] ~ df[[2]], df)
summary(m) %>% print()
lm_eqn_str(df, df[[3]], df[[2]]) %>% cat()
plot1 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[2])) +
geom_boxplot() +
# ggtitle("DV vs IV") +
xlab("IV") + ylab("DV")
plot2 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[3])) +
geom_boxplot() +
# ggtitle("Covar vs IV") +
xlab("IV") + ylab("Covar")
p2 <- ggarrange(plot1, plot2,
labels = c("DV vs IV", "Covar vs IV"),
ncol = 2)
# grid.arrange(plot1, plot2, ncol=2)
print(p2)
cat("\n_____________________\n")
cat("Some stats on the variables\n")
cat("DV - IV:\n")
by(df[[2]], df[[1]], stat.desc) %>% print()
cat("Covar - IV:\n")
by(df[[3]], df[[1]], stat.desc) %>% print()
cat("levene's test\n")
leveneTest(df[[2]], df[[1]], center = median) %>% print()
cat("_____________________\n")
cat("ANOVA of DV with IV\n")
model <- aov(df[[2]] ~ df[[1]], df)
summary(model) %>% print()
cat("_____________________\n")
cat("ANOVA of covariate with IV, independent of DV\n")
model <- aov(df[[3]] ~ df[[1]], df)
summary(model) %>% print()
cat("_____________________\n")
cat("ANOVA of DV with IV, normalized to covariate\n")
c <- df[[2]]/df[[3]]
model <- aov(c ~ df[[1]], df)
summary(model) %>% print()
cat("_____________________\n")
cat("ANCOVA\n")
model <- aov(df[[2]] ~ df[[1]] + df[[3]], data = df)
summary(model) %>% print()
out <- Anova(model, type= "III")
print(out)
return(out)
}
check_ANCOVA <- function(res) {
if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] < 0.05) {
print(paste("Varies significantly with DV and Covar,"), res$`Pr(>F)`[2], res$`Pr(>F)`[3], sep=" ")
} else if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] > 0.05) {
print(paste("Varies significantly with DV,", "p =", res$`Pr(>F)`[2], sep =" "))
} else if (res$`Pr(>F)`[2] > 0.05 && res$`Pr(>F)`[3] < 0.05) {
print(paste("Varies significantly with Covar,", "p =", res$`Pr(>F)`[3], sep=" "))
}
}
Select a tab above to view
total.VO2 <- total %>% select(contains("VO2"))
total.VO2.WT <- total.VO2 %>% select_WT()
total.VO2.KO <- total.VO2 %>% select_KO()
rcf.WT <- lapply(total.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.1))
rcf.KO <- lapply(total.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.1))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in complete cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 1.53037412338406"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 1.64418607579127"
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.1019, num df = 11, denom df = 11, p-value = 0.8751
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3172068 3.8276027
## sample estimates:
## ratio of variances
## 1.101881
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -1.5795, df = 21.948, p-value = 0.1285
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.24455700 0.03310949
## sample estimates:
## mean of x mean of y
## 1.554837 1.660560
df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5032 -0.9561 -0.1400 0.6464 4.4487
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.605 3.295 7.772 9.52e-08 ***
## df[[2]] 1.016 2.039 0.498 0.623
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.654 on 22 degrees of freedom
## Multiple R-squared: 0.01116, Adjusted R-squared: -0.03379
## F-statistic: 0.2482 on 1 and 22 DF, p-value: 0.6233
##
## y = 26 + 1x, r^2 = 0.0112
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.37522815 2.02098097 0.64575282
## sum median mean SE.mean CI.mean.0.95 var
## 19.92672416 1.65187506 1.66056035 0.04616875 0.10161673 0.02557864
## std.dev coef.var
## 0.15993323 0.09631281
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.31808933 1.89266308 0.57457375
## sum median mean SE.mean CI.mean.0.95 var
## 18.65803908 1.56089558 1.55483659 0.04846357 0.10666760 0.02818461
## std.dev coef.var
## 0.16788274 0.10797452
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.7281 0.4027
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0671 0.06707 2.495 0.128
## Residuals 22 0.5914 0.02688
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0002107 2.107e-04 6.008 0.0227 *
## Residuals 22 0.0007716 3.507e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0671 0.06707 2.517 0.128
## df[[3]] 1 0.0319 0.03188 1.196 0.286
## Residuals 21 0.5595 0.02664
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.07817 1 2.9337 0.10147
## df[[1]] 0.09160 1 3.4378 0.07782 .
## df[[3]] 0.03188 1 1.1964 0.28643
## Residuals 0.55952 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
light.VO2 <- light %>% select(contains("VO2"))
light.VO2.WT <- light.VO2 %>% select_WT()
light.VO2.KO <- light.VO2 %>% select_KO()
rcf.WT <- lapply(light.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.1))
rcf.KO <- lapply(light.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.1))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in light cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 1.38025535607356"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 1.44911832907739"
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 0.70246, num df = 11, denom df = 11, p-value = 0.568
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.2022238 2.4401508
## sample estimates:
## ratio of variances
## 0.7024646
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -1.2306, df = 21.348, p-value = 0.2319
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.22122497 0.05663625
## sample estimates:
## mean of x mean of y
## 1.381931 1.464225
df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4395 -1.0950 -0.1280 0.6293 4.5810
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 23.583 2.895 8.145 4.37e-08 ***
## df[[2]] 2.568 2.021 1.271 0.217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.606 on 22 degrees of freedom
## Multiple R-squared: 0.06837, Adjusted R-squared: 0.02602
## F-statistic: 1.614 on 1 and 22 DF, p-value: 0.2171
##
## y = 24 + 2.6x, r^2 = 0.0684
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.19716087 1.91023222 0.71307134
## sum median mean SE.mean CI.mean.0.95 var
## 17.57070207 1.43327793 1.46422517 0.05125155 0.11280390 0.03152066
## std.dev coef.var
## 0.17754058 0.12125224
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.17764983 1.62562865 0.44797882
## sum median mean SE.mean CI.mean.0.95 var
## 16.58316975 1.37143501 1.38193081 0.04295555 0.09454452 0.02214215
## std.dev coef.var
## 0.14880238 0.10767715
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0 0.9945
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0406 0.04063 1.514 0.231
## Residuals 22 0.5903 0.02683
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0001430 1.430e-04 4.812 0.0391 *
## Residuals 22 0.0006536 2.971e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0406 0.04063 1.682 0.2088
## df[[3]] 1 0.0829 0.08287 3.429 0.0782 .
## Residuals 21 0.5074 0.02416
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.01357 1 0.5618 0.46187
## df[[1]] 0.08037 1 3.3260 0.08246 .
## df[[3]] 0.08287 1 3.4294 0.07816 .
## Residuals 0.50743 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
dark.VO2 <- dark %>% select(contains("VO2"))
dark.VO2.WT <- dark.VO2 %>% select_WT()
dark.VO2.KO <- dark.VO2 %>% select_KO()
rcf.WT <- lapply(dark.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.1))
rcf.KO <- lapply(dark.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.1))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in dark cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 1.76003338341844"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 1.8830318193838"
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.4877, num df = 11, denom df = 11, p-value = 0.521
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.4282685 5.1677383
## sample estimates:
## ratio of variances
## 1.487676
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -1.7045, df = 21.186, p-value = 0.1029
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.27797435 0.02748365
## sample estimates:
## mean of x mean of y
## 1.770459 1.895704
df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6653 -0.9232 -0.0703 0.7130 4.0752
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.4603 3.4016 8.367 2.78e-08 ***
## df[[2]] -0.6671 1.8464 -0.361 0.721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.659 on 22 degrees of freedom
## Multiple R-squared: 0.005897, Adjusted R-squared: -0.03929
## F-statistic: 0.1305 on 1 and 22 DF, p-value: 0.7213
##
## y = 28 + -0.67x, r^2 = 0.0059
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.56026454 2.16663334 0.60636880
## sum median mean SE.mean CI.mean.0.95 var
## 22.74845111 1.90175955 1.89570426 0.04658797 0.10253943 0.02604527
## std.dev coef.var
## 0.16138547 0.08513220
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 1.48426845 2.15402819 0.66975974
## sum median mean SE.mean CI.mean.0.95 var
## 21.24550693 1.76912376 1.77045891 0.05682350 0.12506767 0.03874692
## std.dev coef.var
## 0.19684237 0.11118155
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.2641 0.273
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0941 0.09412 2.905 0.102
## Residuals 22 0.7127 0.03240
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0002781 2.781e-04 5.453 0.0291 *
## Residuals 22 0.0011218 5.099e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.0941 0.09412 2.777 0.110
## df[[3]] 1 0.0010 0.00099 0.029 0.866
## Residuals 21 0.7117 0.03389
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.24112 1 7.1143 0.01442 *
## df[[1]] 0.09035 1 2.6657 0.11743
## df[[3]] 0.00099 1 0.0291 0.86615
## Residuals 0.71173 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
Select a tab above to view
total.RQ <- total %>% select(contains("RQ"))
total.RQ.WT <- total.RQ %>% select_WT()
total.RQ.KO <- total.RQ %>% select_KO()
rcf.WT <- lapply(total.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.05))
rcf.KO <- lapply(total.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.05))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in complete cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 0.837304959775152"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 0.847557697979224"
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 2.0112, num df = 11, denom df = 11, p-value = 0.262
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5789754 6.9862557
## sample estimates:
## ratio of variances
## 2.011186
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -1.6606, df = 19.771, p-value = 0.1126
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.035529592 0.004046878
## sample estimates:
## mean of x mean of y
## 0.8336410 0.8493824
df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5901 -0.9979 0.0213 0.8451 4.2127
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.790 12.043 1.726 0.0983 .
## df[[2]] 7.662 14.305 0.536 0.5976
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.653 on 22 degrees of freedom
## Multiple R-squared: 0.01287, Adjusted R-squared: -0.032
## F-statistic: 0.2869 on 1 and 22 DF, p-value: 0.5976
##
## y = 21 + 7.7x, r^2 = 0.0129
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 1.200000e+01 0.000000e+00 0.000000e+00 8.225666e-01 8.818145e-01 5.924789e-02
## sum median mean SE.mean CI.mean.0.95 var
## 1.019259e+01 8.507364e-01 8.493824e-01 5.462712e-03 1.202335e-02 3.580946e-04
## std.dev coef.var
## 1.892339e-02 2.227900e-02
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.000000000 0.000000000 0.000000000 0.787188268 0.870610070 0.083421802
## sum median mean SE.mean CI.mean.0.95 var
## 10.003692068 0.838138435 0.833641006 0.007747015 0.017051066 0.000720195
## std.dev coef.var
## 0.026836449 0.032191853
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.078 0.3104
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.001487 0.0014867 2.758 0.111
## Residuals 22 0.011861 0.0005391
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 1.656e-05 1.656e-05 5.3 0.0312 *
## Residuals 22 6.874e-05 3.124e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.001487 0.0014867 2.804 0.109
## df[[3]] 1 0.000726 0.0007262 1.370 0.255
## Residuals 21 0.011135 0.0005302
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.042923 1 80.9508 1.192e-08 ***
## df[[1]] 0.002041 1 3.8496 0.06315 .
## df[[3]] 0.000726 1 1.3697 0.25498
## Residuals 0.011135 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
light.RQ <- light %>% select(contains("RQ"))
light.RQ.WT <- light.RQ %>% select_WT()
light.RQ.KO <- light.RQ %>% select_KO()
rcf.WT <- lapply(light.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.05))
rcf.KO <- lapply(light.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.05))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in light cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 0.784812573550401"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 0.792340459568411"
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 2.8464, num df = 11, denom df = 11, p-value = 0.09688
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.8194132 9.8875191
## sample estimates:
## ratio of variances
## 2.846395
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -0.56504, df = 17.88, p-value = 0.5791
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.02915937 0.01680362
## sample estimates:
## mean of x mean of y
## 0.7858204 0.7919983
df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3998 -1.1543 -0.2562 1.0993 3.1561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.616 9.587 0.899 0.3786
## df[[2]] 23.604 12.146 1.943 0.0649 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.537 on 22 degrees of freedom
## Multiple R-squared: 0.1465, Adjusted R-squared: 0.1077
## F-statistic: 3.777 on 1 and 22 DF, p-value: 0.06487
##
## y = 8.6 + 24x, r^2 = 0.147
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 1.200000e+01 0.000000e+00 0.000000e+00 7.573150e-01 8.328037e-01 7.548870e-02
## sum median mean SE.mean CI.mean.0.95 var
## 9.503979e+00 7.943118e-01 7.919983e-01 5.574838e-03 1.227013e-02 3.729458e-04
## std.dev coef.var
## 1.931180e-02 2.438364e-02
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.000000000 0.000000000 0.000000000 0.736403245 0.843938825 0.107535580
## sum median mean SE.mean CI.mean.0.95 var
## 9.429844919 0.783804817 0.785820410 0.009405455 0.020701266 0.001061551
## std.dev coef.var
## 0.032581451 0.041461701
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.0908 0.09264 .
## 22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000229 0.0002290 0.319 0.578
## Residuals 22 0.015779 0.0007172
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 9.590e-06 9.594e-06 4.187 0.0529 .
## Residuals 22 5.041e-05 2.291e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.000229 0.000229 0.381 0.544
## df[[3]] 1 0.003169 0.003169 5.277 0.032 *
## Residuals 21 0.012610 0.000600
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.0262812 1 43.7657 1.501e-06 ***
## df[[1]] 0.0010525 1 1.7528 0.19976
## df[[3]] 0.0031690 1 5.2773 0.03199 *
## Residuals 0.0126104 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with Covar, p = 0.03198616938202"
dark.RQ <- dark %>% select(contains("RQ"))
dark.RQ.WT <- dark.RQ %>% select_WT()
dark.RQ.KO <- dark.RQ %>% select_KO()
rcf.WT <- lapply(dark.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.05))
rcf.KO <- lapply(dark.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.05))
data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)
Models <- lapply(data, function(x){
nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})
overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in dark cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))
print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))
## [1] "WT EC50: 0.878325786494703"
print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))
## [1] "KO EC50: 0.903008775451154"
t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])
##
## F test to compare two variances
##
## data: t[, 1] and t[, 2]
## F = 1.2216, num df = 11, denom df = 11, p-value = 0.7458
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3516579 4.2433103
## sample estimates:
## ratio of variances
## 1.221554
t.test(t[,1], t[,2])
##
## Welch Two Sample t-test
##
## data: t[, 1] and t[, 2]
## t = -2.5364, df = 21.783, p-value = 0.01889
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.048071153 -0.004809213
## sample estimates:
## mean of x mean of y
## 0.8730539 0.8994941
df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)
## linear fit of dependent variable with covariates
##
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4826 -1.1494 -0.1335 0.8566 3.8316
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.25 10.47 3.844 0.000882 ***
## df[[2]] -14.68 11.81 -1.243 0.226834
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.608 on 22 degrees of freedom
## Multiple R-squared: 0.06565, Adjusted R-squared: 0.02318
## F-statistic: 1.546 on 1 and 22 DF, p-value: 0.2268
##
## y = 40 + -15x, r^2 = 0.0657
##
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 1.200000e+01 0.000000e+00 0.000000e+00 8.553261e-01 9.413496e-01 8.602341e-02
## sum median mean SE.mean CI.mean.0.95 var
## 1.079393e+01 9.058534e-01 8.994941e-01 6.993827e-03 1.539331e-02 5.869634e-04
## std.dev coef.var
## 2.422733e-02 2.693439e-02
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 1.200000e+01 0.000000e+00 0.000000e+00 8.269545e-01 8.981003e-01 7.114580e-02
## sum median mean SE.mean CI.mean.0.95 var
## 1.047665e+01 8.850517e-01 8.730539e-01 7.729852e-03 1.701329e-02 7.170074e-04
## std.dev coef.var
## 2.677699e-02 3.067049e-02
## Covar - IV:
## df[[1]]: KO
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 28.00000000 3.30000000
## sum median mean SE.mean CI.mean.0.95 var
## 320.70000000 26.80000000 26.72500000 0.28394542 0.62495965 0.96750000
## std.dev coef.var
## 0.98361578 0.03680508
## ------------------------------------------------------------
## df[[1]]: WT
## nbr.val nbr.null nbr.na min max range
## 12.00000000 0.00000000 0.00000000 24.70000000 31.50000000 6.80000000
## sum median mean SE.mean CI.mean.0.95 var
## 333.00000000 27.90000000 27.75000000 0.57689083 1.26972816 3.99363636
## std.dev coef.var
## 1.99840846 0.07201472
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.6755 0.4199
## 22
## _____________________
## ANOVA of DV with IV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.004194 0.004194 6.433 0.0188 *
## Residuals 22 0.014344 0.000652
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 6.30 6.304 2.541 0.125
## Residuals 22 54.57 2.481
## _____________________
## ANOVA of DV with IV, normalized to covariate
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 2.661e-05 2.661e-05 5.489 0.0286 *
## Residuals 22 1.066e-04 4.847e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
## Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]] 1 0.004194 0.004194 6.237 0.0209 *
## df[[3]] 1 0.000220 0.000220 0.327 0.5734
## Residuals 21 0.014124 0.000673
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
##
## Response: df[[2]]
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.068979 1 102.5628 1.552e-09 ***
## df[[1]] 0.003197 1 4.7543 0.04075 *
## df[[3]] 0.000220 1 0.3272 0.57337
## Residuals 0.014124 21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0407456652336419"
reqpkg <- c("ggfortify","rgl","lfda")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "ggfortify"
## [1] '0.4.10'
## [1] "rgl"
## [1] '0.100.54'
## [1] "lfda"
## [1] '1.1.3'
knitr::knit_hooks$set(webgl = hook_webgl)
pca_mutate <- function(data, channel) {
df <- data %>% select(contains(channel))
WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt() %>% select(-variable)
KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt() %>% select(-variable)
df <- bind_cols(WT, KO) %>% set_colnames(c("WT","KO")) %>% melt()
return(df)
}
plot_pca <- function(pca, data, title="") {
autoplot(pca, data = na.omit(data), colour = "id", loadings = TRUE, loadings.colour = "black", loadings.label = TRUE, loadings.label.colour = "black", loadings.label.size = 5) + ggtitle(title)
}
Small filtering step: during the day, there were points where the gas analyzers failed, leaving VO2 and RQ values unbelievably low. The rows contianing data for those points were simply discarded.
light <- light %>% slice(-320:-325)
VO2_light <- pca_mutate(light, "VO2")
RQ_light <- pca_mutate(light, "RQ")
Food_light <- pca_mutate(light, "Food")
Water_light <- pca_mutate(light, "Water")
BM_light <- pca_mutate(light, "BodyMass")
PedMeters_light <- pca_mutate(light, "PedMeters_R")
df_light <- RQ_light %>% set_colnames(c("id","RQ")) %>% add_column(VO2=VO2_light$value) %>% add_column(FoodIn=Food_light$value) %>% add_column(WaterIn=Water_light$value) %>% add_column(BodyMass=BM_light$value) %>% add_column(PedMeters=PedMeters_light$value)
PCA of all variables
df_light %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title = "PCA of all variables, light cycle")
PCA without BodyMass
df_light %>% select(-id,-BodyMass) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA of all variables besides BM, light cycle")
VO2 vs RQ
df_light %>% select(VO2,RQ) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA VO2 and RQ, light cycle")
Food vs Water
df_light %>% select(FoodIn,WaterIn) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA Food and Water, light cycle")
VO2_dark <- pca_mutate(dark, "VO2")
RQ_dark <- pca_mutate(dark, "RQ")
Food_dark <- pca_mutate(dark, "Food")
Water_dark <- pca_mutate(dark, "Water")
BM_dark <- pca_mutate(dark, "BodyMass")
PedMeters_dark <- pca_mutate(dark, "PedMeters_R")
df_dark <- RQ_dark %>% set_colnames(c("id","RQ")) %>% add_column(VO2=VO2_dark$value) %>% add_column(FoodIn=Food_dark$value) %>% add_column(WaterIn=Water_dark$value) %>% add_column(BodyMass=BM_dark$value) %>% add_column(PedMeters=PedMeters_dark$value)
PCA of all variables
df_dark %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, "PCA of all variables, dark cycle")
PCA without BodyMass
df_dark %>% select(-id,-BodyMass) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA of all variables besides BM, dark cycle")
VO2 vs RQ
df_dark %>% select(VO2,RQ) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA VO2 and RQ, dark cycle")
Food vs Water
df_dark %>% select(FoodIn,WaterIn) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA Food and Water, dark cycle")
I attempted to model and plot ldfa. LDFA should give similar results to PCA, with the advantage of being sensitive to multimodality of data, i.e. multiple effects driving an observation, as well as sensitive to distinctiveness in the data that might be lost in reducing dimensions in PCA. In other words, a learning model like LDFA might account for more relationships than a principal component.
More information:
model <- lfda(df_light[-1], df_light[, 1], r=3, metric = "plain")
autoplot(model, data=df_light, frame = T, colour = "id")
plot(model, labels = df_light[, 1])
You must enable Javascript to view this page properly.
model <- lfda(df_dark[-1], df_dark[, 1], r=3, metric = "plain")
autoplot(model, data=df_dark, frame = T, colour = "id")
plot(model, labels = df_dark[, 1])
You must enable Javascript to view this page properly.
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lfda_1.1.3 rgl_0.100.54 ggfortify_0.4.10 WRS2_1.0-0
## [5] pastecs_1.3.21 multcomp_1.4-12 TH.data_1.0-10 MASS_7.3-51.5
## [9] survival_3.1-11 mvtnorm_1.1-0 effects_4.1-4 compute.es_0.2-4
## [13] car_3.0-7 carData_3.0-3 nplr_0.1-7 purrr_0.3.3
## [17] scales_1.1.0 tibble_2.1.3 docstring_1.0.0 reshape2_1.4.3
## [21] DT_0.13 ggpubr_0.2.5 ggplot2_3.3.0 dplyr_0.8.5
## [25] magrittr_1.5
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_1.4-1 ggsignif_0.6.0
## [4] rio_0.5.16 htmlTable_1.13.3 markdown_1.1
## [7] base64enc_0.1-3 mc2d_0.1-18 rstudioapi_0.11
## [10] roxygen2_7.1.0 farver_2.0.3 RSpectra_0.16-0
## [13] xml2_1.2.5 codetools_0.2-16 splines_3.6.3
## [16] knitr_1.28 Formula_1.2-3 jsonlite_1.6.1
## [19] nloptr_1.2.2.1 cluster_2.1.0 png_0.1-7
## [22] shiny_1.4.0.2 compiler_3.6.3 backports_1.1.5
## [25] fastmap_1.0.1 assertthat_0.2.1 Matrix_1.2-18
## [28] survey_3.37 later_1.0.0 acepack_1.4.1
## [31] htmltools_0.4.0 tools_3.6.3 gtable_0.3.0
## [34] glue_1.3.2 Rcpp_1.0.4.6 cellranger_1.1.0
## [37] vctrs_0.2.4 nlme_3.1-144 crosstalk_1.1.0.1
## [40] xfun_0.12 stringr_1.4.0 openxlsx_4.1.4
## [43] lme4_1.1-21 miniUI_0.1.1.1 mime_0.9
## [46] lifecycle_0.2.0 klippy_0.0.0.9500 zoo_1.8-7
## [49] hms_0.5.3 promises_1.1.0 sandwich_2.5-1
## [52] RColorBrewer_1.1-2 yaml_2.2.1 curl_4.3
## [55] gridExtra_2.3 rpart_4.1-15 reshape_0.8.8
## [58] latticeExtra_0.6-29 stringi_1.4.6 checkmate_2.0.0
## [61] manipulateWidget_0.10.1 boot_1.3-24 zip_2.0.4
## [64] rlang_0.4.5 pkgconfig_2.0.3 evaluate_0.14
## [67] lattice_0.20-40 htmlwidgets_1.5.1 labeling_0.3
## [70] cowplot_1.0.0 tidyselect_1.0.0 plyr_1.8.6
## [73] R6_2.4.1 Hmisc_4.4-0 DBI_1.1.0
## [76] pillar_1.4.3 haven_2.2.0 foreign_0.8-75
## [79] withr_2.1.2 mgcv_1.8-31 abind_1.4-5
## [82] nnet_7.3-13 crayon_1.3.4 rARPACK_0.11-0
## [85] rmarkdown_2.1 jpeg_0.1-8.1 grid_3.6.3
## [88] readxl_1.3.1 data.table_1.12.8 forcats_0.5.0
## [91] webshot_0.5.2 digest_0.6.25 xtable_1.8-4
## [94] tidyr_1.0.2 httpuv_1.5.2 munsell_0.5.0
## [97] mitools_2.4